/****************************************************************************
 *                                                                          *
 * File    : swhmodel.c                                                     *
 *                                                                          *
 * Purpose : CPRE588 Solar Water Heater Model.                              *
 *                                                                          *
 * History : Date      Reason                                               *
 *           04/18/11  Created                                              *
 *                                                                          *
 ****************************************************************************/
#include <math.h>
#include <string.h>
#include <stdio.h>

#include "solpos00.h" 

   struct posdata pd, *pdat; /* declare a posdata struct and a pointer for
                                 it (if desired, the structure could be
                                 allocated dynamically with malloc) */
   long retval;              /* to capture S_solpos return codes */
 
   long  secondsSinceMidnight ;
   long  checkInterval ;
   long  currentTime ;
   int   s_sec, s_min, s_hour, s_temp;
   long  startTime, finishTime ;
   
   float sunAzimuth, sunElevation;
   float colAzimuth, colElevation;
   float azimuthStepSize, elevationStepSize;
   float diodeOffset;
   float pointingError;
   
   int   deltaTemp;
   int   desiredWaterTemp;
   int   lastWaterTemp;
   int   lastTankTemp ;
   int   actualWaterTemp;
   int   colWaterTemp;
   
   enum {pumpOFF, pumpSLOW, pumpMEDIUM, pumpFAST} pumpSpeed; 

void setCollectorAzimuth(float desiredAzimuth);
void setCollectorElevation(float desiredElevation);
float getDiodeResponse(float incidentAngle);
void setPumpSpeed(void);
   
/****************************************************************************
 *                                                                          *
 * Function: main                                                           *
 *                                                                          *
 * Purpose : Main function.                                                 *
 *                                                                          *
 ****************************************************************************/
int main(int argc, char *argv[])
{
   /* Initialize system variables */
   /*-----------------------------*/   
   checkInterval = 60*5; // Check every 5 minutes

   /* Initialize temperatures */
   actualWaterTemp  = lastWaterTemp = 75  ; /* Tank water temperature */
   lastTankTemp     = actualWaterTemp ;
   desiredWaterTemp = 140 ;   /* Desired temp of hot water in tank */ 
   colWaterTemp     = 75  ;   /* Initial temp of water in dish collector */

   /* Resolution of stepper motors, in degrees */
   azimuthStepSize   = 1.2 ;
   elevationStepSize = 1.2 ;

   /* Angle offset of photo diodes from boresight of dish */
   diodeOffset = 45.0 ;

   /* Pump is initially turned off */
   pumpSpeed = pumpOFF;
    
   pdat = &pd; /* point to the structure for convenience */

    /* Initialize structure to default values. (Optional only if ALL input
       parameters are initialized in the calling code, which they are not
       in this example.) */
   S_init(pdat);

   pdat->function &= ~S_DOY ; /* (sets month and day input) */

    /* Use Fredericksburg, VA */
   pdat->longitude = -77.457283;  /* Note that latitude and longitude are  */
   pdat->latitude  =  38.356778;  /*   in DECIMAL DEGREES, not Deg/Min/Sec */
   pdat->timezone  =  -5.0;   /* Eastern time zone, even though longitude would
                                  suggest Central.  We use what they use.
                                  DO NOT ADJUST FOR DAYLIGHT SAVINGS TIME. */

   pdat->year     = 2011;   /* The year is 2011. */
   pdat->month    =   5;    /* Set month to May. */
   pdat->day      =   6;    /* May 6th, deadline for this project.*/

   /* The time of day (STANDARD time) is 4 AM */
   pdat->hour   =  4 ;
   pdat->minute =  0 ;
   pdat->second =  0 ;

   retval = S_solpos (pdat);  /* S_solpos function call */
   S_decode(retval, pdat);    /* ALWAYS look at the return code! */
  
  /* Print initial conditions   */  
   printf("Solar Water Heater\n");
   printf("==================\n");
   printf(" Desired Temp : %d\n",desiredWaterTemp);
   printf(" Starting Temp: %d\n",actualWaterTemp);
   printf(" -----------------\n");       
   printf(" Latitude : %f\n",pdat->latitude);
   printf(" Longitude: %f\n",pdat->longitude);   
   printf(" Date: %d.%0.2d.%0.2d\n",pdat->year, pdat->month, pdat->day);

   /* Determine the sunrise time  */
   /*-----------------------------*/  
   secondsSinceMidnight = (long) (pdat->sretr * 60);
   s_hour = (int) (secondsSinceMidnight / 3600);   /* hours since midnight [0,23] */
   s_temp = s_hour*3600;
   s_min  = (int) ( secondsSinceMidnight - s_temp );
   s_min  = (int) (s_min / 60);                    /* minutes after the hour [0,59] */        
   s_sec  = (int) (secondsSinceMidnight % 60);     /* seconds after the minute [0,60] */
   printf ( "Sunrise is at: %02d:%02d:%02d\n", s_hour, s_min, s_sec );

   /* Round seconds */
   if( s_sec > 29 ) ++s_min;
   s_sec=0;
   startTime  = (s_hour*3600) + (s_min*60) ; // in seconds   

   /* Set up variables for subsequent function call for initial  */
   pdat->hour   =  s_hour ;
   pdat->minute =  s_min ;
   pdat->second =  s_sec ;

   /* Determine the sunset time   */
   /*-----------------------------*/     
   secondsSinceMidnight = (long) (pdat->ssetr * 60);
   s_hour = (int) (secondsSinceMidnight / 3600);   /* hours since midnight [0,23] */
   s_temp = s_hour*3600;
   s_min  = (int) ( secondsSinceMidnight - s_temp );
   s_min  = (int) (s_min / 60);                    /* minutes after the hour [0,59] */        
   s_sec  = (int) (secondsSinceMidnight % 60);     /* seconds after the minute [0,60] */
   printf ( "Sunset  is at: %02d:%02d:%02d\n", s_hour, s_min, s_sec );

   /* Round seconds */
   if( s_sec > 29 ) ++s_min;
   s_sec=0;
   finishTime  = (s_hour*3600) + (s_min*60) ; // in seconds 

   /* Initialize the collector pointing position */
   retval = S_solpos(pdat);   /* S_solpos function call */
   S_decode(retval, pdat);    /* ALWAYS look at the return code! */
   colAzimuth   = pdat->azim ;
   colElevation = pdat->elevref ;

   /* Print initial conditions   */
   printf("          |         Sun         |        Dish         |      | Tank |\n");   
   printf("   Time   | Azimuth | Elevation | Azimuth | Elevation | Pump | Temp |\n");
   printf("----------|---------|-----------|---------|-----------|------|------|\n");

   /* Start the check loop        */
   /*-----------------------------*/ 
   for(currentTime=startTime;currentTime<finishTime;currentTime+=checkInterval)
   {
    /* For each check interval do the following: */
    /* 1. Calculate the current position of the sun. */
    /* 2. Calculate the tank water temperature based on the last pump speed setting  */
    /*     and the last measured temperature of the water in the tank.               */
    /* 3. Account for usage by reducing hot water tank temperature at certain times  */
    /*     during the day.                                                           */
    /* 4. Compare desired tank water temperature versus actual water temperature.    */
    /*     If the tank water temperature is less than desired move the dish to       */
    /*       point at the sun (heat water in collector system) and turn on the pump  */
    /*       to circulate water into hot water tank heat exchanger, thus heating the */
    /*       tank water.                                                             */
    /*     If the tank water temperature is higher than the desired temperature,     */
    /*       Turn off the pump and lock the dish in its current position (to keep    */
    /*       the water in the collector system from heating.)                        */
    /*===============================================================================*/
    
    /*-----------------------------------------------*/
    /* 1. Calculate the current position of the sun. */
    /*-----------------------------------------------*/
    
    /* Convert the check time into format used by solpos */
      s_hour = (int) ( currentTime / 3600);   /* hours since midnight [0,23] */
      s_temp = s_hour*3600;
      s_min  = (int) ( currentTime - s_temp );
      s_min  = (int) ( s_min / 60);           /* minutes after the hour [0,59] */        
      s_sec  = (int) ( currentTime % 60);     /* seconds after the minute [0,60] */
      
      pdat->hour    =  s_hour;
      pdat->minute  =  s_min;
      pdat->second  =  s_sec;

      /* Calculate the current position of the sun */
      retval = S_solpos (pdat);  /* S_solpos function call */
      S_decode(retval, pdat);    /* ALWAYS look at the return code! */
      
      sunAzimuth   = pdat->azim    ;
      sunElevation = pdat->elevref ;
   
      printf("%02d:%02d:%02d   %6.2f    %6.2f",
         s_hour, s_min, s_sec, sunAzimuth, sunElevation );

      /*----------------------------------------------------------------------------*/
      /* 2. Calculate tank water temperature based on past pump speed and interval. */
      /*----------------------------------------------------------------------------*/

      switch( pumpSpeed )
      {
         case pumpOFF:
            /* Assume exponential decay of tank water temperature */
            /* Actual decay function would be based on efficiency of tank insulation */
            actualWaterTemp = lastWaterTemp*exp(-0.004*(checkInterval/60));
            /* Never let tank water temperature get less than ambient */
            /* Assumed to be 75 degrees F */
            if( actualWaterTemp < 75 ) actualWaterTemp = 75;
            break;

        /* Very simple model for heat transfer */
        /* Assume a linear transfer of heat from collector water to tank water. */
        /* Faster pump speed equates to a faster transfer of heat */
         case pumpSLOW:
            actualWaterTemp = lastWaterTemp + (1*(checkInterval/60));
         break;
         
         case pumpMEDIUM:
            actualWaterTemp = lastWaterTemp + (3*(checkInterval/60));
         break;
         
         case pumpFAST:
            actualWaterTemp = lastWaterTemp + (6*(checkInterval/60));
         break;
         
         default:

         break;
      }
      lastWaterTemp = actualWaterTemp;

      /*------------------------------------------------*/
      /* 3. Model hot water usage during daytime hours. */
      /*------------------------------------------------*/
      /* Simple assumption is tank water temp returns to ambient  */
      /* Account for usage at certain times of day                */
      /* 9:30AM  = 9*3600 + 30* 60 = 34200 seconds after midnight */
      /* 12:00PM = 12*3600 =  43200 seconds after midnight        */
      /*  3:30PM = 15*3600 + 30*60 = 55800 seconds after midnight */
      
      if( (currentTime >= 33000) && (33500 >= currentTime ) )
      {
         actualWaterTemp = lastWaterTemp = 75;
      }
      else if( (currentTime >= 42000) && (42500 >= currentTime) )
      {
         actualWaterTemp = lastWaterTemp = 75;

      }
      else if( (currentTime >= 54000) && (54500 >= currentTime) )
      {
         actualWaterTemp = lastWaterTemp = 75;
      }

      /*-------------------------------------------------------*/
      /* 4. Compare tank water temp against user desired temp. */
      /*-------------------------------------------------------*/
      
      /* Check hot water tank temperature versus desired temperature */
      if( actualWaterTemp < desiredWaterTemp )
      { /* Need to heat water in tank */
         /* Need to heat the collector water, point dish at sun */
         setCollectorAzimuth(pdat->azim);
         setCollectorElevation(pdat->elevref);
         
         /* Calculate water temp of solar collector water based from ave pointing error */
         /* Assumes a simple linear heat transfer function based on angle of dish */
         pointingError = fabs( ( (pdat->azim - colAzimuth)/pdat->azim)*100.0);
         pointingError += (fabs( ( (pdat->elevref - colElevation)/pdat->elevref) ) *100.0);
         pointingError /= 2 ;
         colWaterTemp = (int)(-1.37*pointingError + 212);

         /* Set pump speed based on temperature differential from desired temp and actual temp */
         if( actualWaterTemp > desiredWaterTemp )
         {
            pumpSpeed = pumpOFF;
         }
         else
         {
            deltaTemp = desiredWaterTemp - actualWaterTemp;
            if( deltaTemp > 35 )
            {
               pumpSpeed = pumpFAST;
            }
            else if ( deltaTemp > 10 )
            {
               pumpSpeed = pumpMEDIUM;
            }
            else
            {
               pumpSpeed = pumpSLOW;
            }
         }
      }
      else
      { /* Tank water temperature is higher than desired */
        /* Dish position is not updated; i.e. it's locked in place */
        /* Turn off the pump */
         pumpSpeed = pumpOFF;
      }
      
      printf("       %6.2f   %6.2f",colAzimuth, colElevation);

      switch( pumpSpeed ){

         case pumpOFF:  printf("        OFF");
         break;
         
         case pumpSLOW: printf("       SLOW");
         break;
         
         case pumpMEDIUM:printf("        MED");
         break;
         
         case pumpFAST:printf("       FAST");
         break;
         
         default:
         printf("        ERR");
         break;
      }
      printf("    %d\n",actualWaterTemp);
         
   }

  printf(" DONE\n");
  return 0;
}

/****************************************************************************
 *                                                                          *
 * Function: setCollectorAzimuth                                            *
 *                                                                          *
 * Purpose : Sets the azimuth direction of the solar collector, as close as *
 *           possible, based on the feedback from photodiode detectors and  *
 *           the angular resolution of the stepper motor.                   *
 *                                                                          *
 ****************************************************************************/
void setCollectorAzimuth(float desiredAzimuth)
{
   float deltaAngle       ;
   float diodeX1          ;
   float diodeX2          ;
static float collectorAzimuth ;
static float deadband  = 0.1 ;
      
   /* What is current azimuth setting of stepper motor? */
   collectorAzimuth = colAzimuth;

   /* Initialize the delta angle */
   deltaAngle = desiredAzimuth - collectorAzimuth;

   /* Determine the new setting for the azimuth stepper motor */
   while( fabs(deltaAngle) > azimuthStepSize )
   {
      /* Fetch the voltage response of each photodiode */         
      diodeX1 = getDiodeResponse(diodeOffset+deltaAngle);
      diodeX2 = getDiodeResponse(diodeOffset-deltaAngle);       

      /* Move the dish based on the response from the diodes */
      /* Assumes setting the stepper motor is as simple as writing a value to */
      /* a memory location.                                                   */ 
      if( (diodeX1 > diodeX2) && ((diodeX1 - diodeX2) > deadband) )
      {
         collectorAzimuth -= azimuthStepSize; /* Move dish left */
      }
      else if( (diodeX2 > diodeX1) && ((diodeX2 - diodeX1) > deadband) )
      {
         collectorAzimuth += azimuthStepSize; /* Move dish right */
      }

      /* Determine how close to desired azimuth */
      deltaAngle = desiredAzimuth - collectorAzimuth;
   }
  
  /* Save azimuth setting to a global variable */
  colAzimuth = collectorAzimuth;
}

/****************************************************************************
 *                                                                          *
 * Function: setCollectorElevation                                          *
 *                                                                          *
 * Purpose : Sets the elevation of the solar collector, as close as         *
 *           possible, based on the feedback from photodiode detectors and  *
 *           the angular resolution of the stepper motor.                   *
 *                                                                          *
 ****************************************************************************/
void setCollectorElevation(float desiredElevation)
{
   float deltaAngle      ;
   float diodeYpos       ;
   float diodeYneg       ;
static float collectorElevation ;
static float deadband  = 0.1 ;
 
  /* What is current elevation setting of stepper motor? */
  collectorElevation = colElevation;

  /* Initialize the delta angle */
  deltaAngle = desiredElevation - collectorElevation;

   /* Determine the new setting for the elevation stepper motor */
   while( fabs(deltaAngle) > elevationStepSize )
   {
      /* Fetch the voltage response of each photodiode */ 
      diodeYpos = getDiodeResponse(diodeOffset-deltaAngle);
      diodeYneg = getDiodeResponse(diodeOffset+deltaAngle);       

      /* Move the dish based on the response from the diodes */
      /* Assumes setting the stepper motor is as simple as writing a value to */
      /* a memory location.                                                   */ 
      if( (diodeYpos > diodeYneg) && ((diodeYpos - diodeYneg) > deadband) )  
      {
         collectorElevation += elevationStepSize; /* Move dish up */
      }
      else if( (diodeYneg > diodeYpos) && ((diodeYneg - diodeYpos) > deadband) )
      {
         collectorElevation -= elevationStepSize; /* Move dish down */
      }
      
      /* Determine how close to desired elevation */
      deltaAngle = desiredElevation - collectorElevation;
   }
   
  /* Save elevation setting to a global variable */
  colElevation = collectorElevation;
}

/****************************************************************************
 *                                                                          *
 * Function: getDiodeResponse                                               *
 *                                                                          *
 * Purpose : Returns the voltage from a photodiode based on the angle of    *
 *           sunlight incident on the detector surface. Return value        *
 *           represents a voltage between 0 and 3.5 volts. The response of  *
 *           the diode is based on a Gaussian function that attempts to     *
 *           model the response of a Texas Instruments TSL25x series        *
 *           light-to-voltage sensor.                                       * 
 *                                                                          *
 ****************************************************************************/
float getDiodeResponse(float incidentAngle)
{
   float temp;
   temp = -(incidentAngle*incidentAngle)/2187.0;
   return (3.5*exp(temp));
}

